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Abstract 

We obtain new results on the stability of discrete dark solitons bifurcating from the anti- 
continuum limit of the discrete nonlinear Schrodinger equation, following the analysis of our 
previous paper [Physica D 212, 1-19 (2005)]. We derive a criterion for stability or instability of 
dark solitons from the limiting configuration of the discrete dark soliton and confirm this criterion 
numerically. We also develop detailed calculations of the relevant eigenvalues for a number of 
prototypical configurations and obtain very good agreement of asymptotic predictions with the 
numerical data. 

In this paper, we address the dynamical lattice model governed by the discrete nonlinear Schrodinger 
(DNLS) equation [5]. We consider the defocusing version of this equation in the form 

iu n + e (n n+ i - 2u n + u n _i) - \u n \ 2 u n = 0, (1) 

where n G Z, u n (t) : R — > C, and e > 0. The stationary solutions u n (t) = 4> n e~ lt are found from 
second-order difference equation 

(<g -l)0 n = e(0 n+ i-20 n + n _i) (2) 

for a real- valued sequence {4> n }nez, denoted in vector notations by (f>. Our aim here is to study 
discrete dark solitons which are defined by the non-vanishing boundary conditions at infinity, e.g. 
lim n _± OQ <j) n = ±1. The limiting configuration of dark solitons at e = is defined by the decom- 
position Z = U + U U- U Uq such that <f> n = ±1 for n G U± and (j) n = for n G Uq. Our previous 
work [10] addressed stability of discrete bright solitons when dim(C/+ U U-) < oo. In this paper, 
we shall study stability of discrete dark solitons, when dim(C/o) < oo and there exists N > 1 such 
that n G U± for all ±n > N. These solutions were considered recently in [31 [12], as well as earlier 
in [U [6] using predominantly numerical computations. 

The topic of dark solitons and their stability is not only of theoretical and mathematical interest, 
but is also a subject of relevance to presently available experimental settings. In particular, discrete 
dark solitons have been observed in the context of AlGaAs waveguide arrays in the anomalous 
diffraction regime [8]. Furthermore, as was illustrated in [3J, similar phenomenology can be observed 
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in the discrete dark solitons that arise in defocusing lithium niobate waveguide arrays which exhibit 
a saturable nonlinearity due to the photovoltaic effect; in the latter case, experimental results were 
presented in the work of [TT]. Although these nonlinear optics experiments are the most promising 
realizations of discrete dark solitons, such waveforms may also be relevant to the atomic physics. In 
particular, dark solitons were considered for Bose-Einstein condensates in the presence of a periodic, 
so-called optical lattice, potential [H [9] (although in the latter setting, discrete dark solitons have 
not yet been experimentally realized). 

To determine the persistence and stability of discrete dark solitons, we study spectra of the linear 
operators L + and L_ defined by 

(L + 1p) n = (3(j)l - l)tp n - € (lpn+1 - 2lpn + tpn-l) , 
(L_V>)n = (4>l, ~ l)lp n ~ £ (lpn+1 ~ 2 VVi + lpn-l) ■ 

If G Z°°(Z) for any e > 0, then the operators L± map / 2 (Z) to itself. Their spectrum at e = is 
computed explicitly. The operator L + has an eigenvalue 2 of multiplicity dim(£/+) +dim(Z7_) = oo 
and an eigenvalue —1 of multiplicity dim(f7o) < oo. The operator L_ has an eigenvalue of 
multiplicity dim(£7 + ) + dim(E/_) = oo and an eigenvalue —1 of multiplicity dim(C/o) < cxd. 

Since l 2 {1) is a Banach algebra with respect to the pointwise multiplication and the operator 
L + is continuously invertible in / 2 (Z) for sufficiently small e > 0, persistence of solutions of the 
difference equation in Z 2 (Z) C /°°(Z) with respect to small parameter e is proved using the 
Implicit Function Theorem. Analysis of the stability problem 

(L + u) n = -Xw n , (L-w) n = An n (3) 

for small e > is, however, more complicated because of the splitting of the zero eigenvalue of 
infinite multiplicity into a spectral band located at 

A s = | A G C : -2i/2e(l + 2e) < ImA < 2y/2e(l + 2e)| 

and a number of isolated eigenvalues around the origin. We shall count these eigenvalues by using 
the recent results of [21 [10] . 

Since (L_0) n = and the non-decaying sequence {</> n }nez does not oscillate as n — ► ±oo, 
is at the bottom of the continuous spectrum of L_, which is located for A G [0, 4e]. By the 
discrete Sturm theory [7], the number of negative eigenvalues of L_ equals the number of times 
the sequence {^ n }nez changes sign on Z. To compute this number, we subdivide Uo into segments 
Uq = UjL^nJ, rij] for some N < oo and denote the number of sign-changes between adjacent nodes 
in U + U U- by gq. 

Lemma 1 There exists eo > such that, for any e G (0, eo), the number of sign changes of 
the sequence {4> n }n£Z equals dim(Z7o) + o~o + Y^j=\ a j> where o~j is associated with the segment 
Uj = [nj,rij] C Uq, such that 

1 if dim([/j) is odd and sign(</> n -_ 1 4> n + + i) = 1 

j j 

if dim(Uj) is odd and sign(</) n -_ 1 </) n + +1 ) = —1 
°^ 1 if dim(Uj) is even and sign((/> n -_ 1 n + +1 ) = —1 ^ 

if dim(C/j) is even and sign(0 n - _ 1 n + +1 ) = 1 
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Proof. By persistence of solutions in Z°°(Z)-norm for sufficiently small e, any sign change between 
the adjacent nodes in U+ U U- persists in e. Therefore, the statement of the lemma is proved if we 
can prove for a particular segment Uj = [rij , rij] of length rij = n~j" — nj + 1 that the number of 
sign changes equals nj + o~j, where Oj is given by To do this with an application of Lemma 2.3 
from [10], we use the staggering transformation <f> n = (—l) n ip n and rewrite the different equation 
(|2|) in the form 

(1 - ^D'Pn = 6 + ty n + </?n-l) • 

n + +l 

By Lemma 2.3 of [10], there is only one sign difference in the sequence {<p n } i if sign(</? +,,) = 

— 1 and none if sign(w - + , J) = 1. If rij is odd, the staggering transformation gives (nj + 1) 

j j 

n + +l 

sign differences in the sequence {(fi n } - if sign(</> n -_i<f> „+ , i) = 1 and n, sign differences if 

ra + +l 

sign(^) n -_-\<t> n + 1 1) = —1- If % is even, however, the sequence {4> n } - has rij sign differences 
if sign(0 -_ 1 <p + , i) = 1 and (rij + 1) sign differences if sign(</> n -_ 1 <p n + + i) = — 1- Thus, formula 
(UJ) is proved. ■ 

Corollary 2 T/ie number Nq = ao + ]Cj=i 17 j equals the number of small negative eigenvalues of 
L_ for e > bifurcating from the zero eigenvalue of infinite multiplicity for e = 0. 

Proof. This follows from the fact that L_ has dim(?7o) negative eigenvalues at e = 0. ■ 

Definition 3 The stability problem (dp is said to have a purely imaginary eigenvalue of negative 
Krein signature i/(L_u, u) = (L^ 1 w, w) < for the corresponding eigenvector (u, w). 

Theorem 4 There exists eo > such that, for any e E (0, eo), the stability problem has exactly 
dim(C/o) pairs of purely imaginary isolated eigenvalues of negative Krein signature bounded away 
from the continuous spectrum and exactly Nq pairs of small real eigenvalues. 

Proof. Since L + is invertible for sufficiently small e, we rewrite the stability problem (|3|) in the 
form 

L_w = 7L^ 1 w, 7 = —A 2 . (5) 
Since L_ is not invertible for any e > 0, we shift the eigenvalue problem to the form 

(L_ + 5L- 1 ) w = ( 7 + 5)L^w, 

for sufficiently small 5 > 0. Since Null(L_) lies in the positive subspace of L+ 1 at e = and the 
number of negative eigenvalues of L_ is unchanged in e G (0, eo), for a fixed e G (0, eo), there is 
a small 5 = 5(e), such that the number of negative eigenvalues of L_ + SL^ 1 is the same as that 
of L_ . Conditions of Theorem 1 of [2] are now satisfied and we count the negative eigenvalues of 
L- + 5L^_ (same as for L_) and (same as for L+) as follows: 

n(L_ ) = N- + N+ + iV c+ , n(L + ) = N~ + N+ + N c+ , 
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where n(L±) denotes the number of negative eigenvalues of L±, N c + denotes the number of complex 
eigenvalues 7 in the upper half-plane, denote the number of positive/negative eigenvalues 7 
with (w, L^ 1 w) < for corresponding eigenvectors w, and N~ denotes the number of negative 
eigenvalues 7 with (w,L^ 1 w) > for corresponding eigenvectors w. By Lemma Q] and Corollary 
[21 we have n(L_) = dim(JJo) + No and n(L + ) = dim(£/o) for sufficiently small e G (0, eo). 

At e = 0, there exists dim([/o) eigenvalues 7 = 1 with (w, i/j^w) < 0, where the sequence 
{wn}nez for the eigenvector w is compactly supported in Uq. By Proposition 5.1 in [2j, the 
eigenvalue 7 = 1 is hence semi-simple (that is algebraic and geometric multiplicities coincide). 
Therefore, all dim(tTo) eigenvalues persist for positive values of 7 for sufficiently small e. By 
continuity of the eigenvectors w in e, the inequality (w, w) < holds for sufficiently small 
e > 0, and therefore, n(L + ) = dim(C/o) = N£, such that N c + = N~ = and N~ = Nq. Therefore, 
all No bifurcations of small negative eigenvalues of L_ for e G (0, eo) from the zero eigenvalue of 
L_ for e = result in pairs of small real eigenvalues A = ±^—7 of the stability problem ([3]). ■ 

Remark 5 Small eigenvalues of the operator L_ can be found from the difference eigenvalue prob- 
lem 

00 

V n ^n - e (^ n+1 + Vn-1 " ^n) = V n = V® + ^ ^ V n ] » (6) 

k=l 

where K (0) = {^nf - 1, VP = 2</4° ) tf ) , K (2) = 2^ O) 0n } + (^V and so on, due to analytic 
dependence of the solution cj> of the difference equation p|) on e. // w is supported in U + U {/_ 
and e = 0, i/ien L+w = 2w. Since I 2 -eigenvectors of the difference equation (0|) /or small negative 
eigenvalues [i are supported in U+ U U- as e —* 0, a small negative eigenvalue fx for L_ is related 
to a small negative eigenvalue 7 /or ff/iaf is £/ie eigenvalue of the stability problem (0) u>ii/i 

f/ie l 2 -eigenvector) by the asymptotic approximation lim^o 7/V = 2. 

As the simplest application of our results, we consider two basic configurations of discrete dark 
solitons from [3j HI [6] . 

• If U± = 7L± and Uq = {0} (a so-called on-site dark soliton), then Nq = and, according to 
Theorem 01 the dark soliton is stable with a single pair of purely imaginary eigenvalues of 
negative Krein signature near A = ±i. 

• If U+ = Z+, U- = Z_ U {0}, and Uq = (a so-called inter-site dark soliton), then Nq = 
ctq = 1 and the dark soliton is unstable with a single pair of real eigenvalues. The asymptotic 
approximation of the unstable eigenvalue can be obtained with the following argument. The 
solution of the difference equation ([2]) is expanded in the power series cj> = 4>^ +^(j> +0(e 2 ), 
where (j)^ is compactly supported with cft^ = 1, 0^ = —1 and (f>n = for all n £ Z\{0, 1}. 
Since Vn°^ = for all n 6 Z and = —2 for n = {0, 1} and Vn = otherwise, the potential 
V of the discrete Schrodinger equation © is negative at the leading order. By the discrete 
Sturm theory [7j, it traps a unique negative eigenvalue with the symmetric eigenfunction 
ip n = ip-n+i, n £ N. Using the parametrization 

fx = e (2 - e K - e~ K ) (7) 
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Figure 1: Numerical approximations of the unstable eigenvalue for the inter-site dark soliton (solid 
line) together with the asymptotic prediction A = \J~^ (dashed curve) and the approximation 
A = y/2e of [3] (thick dashed curve). 



and solving the eigenvalue problem for the eigenvector fa = 1, fa = Ce~ K ( n ~ 2 ^ for n > 2, 
we obtain C = e~ K and e K = 3 at the leading order of O(e), which gives /i = — |e + 0(e 2 ). 
Using Remark [5j we conclude that the pair of real eigenvalues of the stability problem ([3]) 

is given by A = ±y/—j = ±-i/y(l + 0(e)). This approximation is shown on Fig. [1] with 
thin dashed line, while the solid line shows results of numerical approximations of eigenvalues 
of the truncated linear stability system ([3]). It should be noted that the earlier work of [3] 
approximated the real eigenvalue pair of the inter-site dark soliton as A = ±\/2e- As can be 
readily observed from a solid dashed line in Fig. [TJ this asymptotic prediction is incorrect. 



To illustrate more complicated applications of Theorem 01 we consider several composite discrete 
dark solitons, some of which were studied in |12j . 



If U+ = {0} U Z+\{1}, U = 0, and U- = Z_ U {1}, then N = a = 3, such that three 
pairs of real (unstable) eigenvalues occurs in the linearized problem (|3|). To find asymptotic 
approximations of these eigenvalues, we again consider eigenvalues of L_ in the difference 
equation <^ with the potentials = for all n € Z and 



We construct two symmetric eigenvectors and one anti-symmetric eigenvector for three neg- 
ative eigenvalues \i. For symmetric eigenvectors, ip n = V'-n+i; n £ Z with t/ji = 1, fa = B, 
Tp n = Ce~ K( - n - 3 \ n > 3, we use the parametrization ([7]) and obtain at the leading order of 
O(e): 

C = Be- K , (e K -2)B = l, B = e K + e~ K - 5. 



1, 


n = —1 








-2, 


n = 




f " 4 ' 


n = {0,l} 


2, 


n = 1 , 


= 20(M 1} = < 


-2, 


n = {-l,2}, 


-1, 


n = 2 




I 0, 


otherwise 


0, 


otherwise 
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Figure 2: The left panel compares the asymptotic predictions (dashed lines) with the results of 
numerical linear stability analysis (solid lines) for the three positive real eigenvalues of the discrete 
dark soliton with U + = {0} U Z + \{1}, U$ = 0, and £/_ = Z_ U {1}. The right panel shows a 
typical example of the solution profile (top) and the corresponding spectral plane A = A r + i\ of 
its linearization spectrum for e = 0.05. 



Eliminating B, we obtain a cubic equation for z = e K : 

z 3 - 7z 2 + lOz - 2 = 0. (8) 

There exist two solutions of the cubic equation in the interval z > 1, namely z\ ~ 1.63667 
and zi ~ 5.12489. For the anti-symmetric eigenvector, ifj n = — tp_ n+ i, n £ Z with tpi = 1, 
ip2 = B, ip n = Ce~ K ^ n ~ 3 \ n > 3, we obtain at the leading order of 0(e): 

C = Be~ K , {e K -2)B = l, B = e K + e~ K -3. 

Eliminating B, we obtain a cubic equation for z = e K : 

z 3 - 5z 2 + 6z - 2 = 0. (9) 

Since one root is z = 1, we can find a unique root of the cubic equation in the interval 
z > 1, namely z% = 2 + y/2. Each of the three roots above generates a negative eigenvalue 
for fx = e(2 — z — z^ 1 ) + 0(e 2 ). Each negative eigenvalue ^ of L_ generates a negative 
eigenvalue 7 of L + L_ with the correspondence 7 = 2fi + 0(e 2 ). Summarizing, the three 
pairs of real eigenvalues are given asymptotically by A = ±0.70380-^6, A = ±1.84776- v /e 
and A = ±2.57683y / e. These theoretical predictions are compared with the results of full 
numerical linear stability analysis in the left panel of Fig. [21 yielding a good agreement for 
small values of e. A typical example of the discrete dark soliton and its linearization spectrum 
for e = 0.05 is shown in the right panel of Fig. [2j 

• If U + = Z+\{1}, U Q = {0}, and U- = Z_ U {1}, then a = 1, 01 = 1, such that = 2 and 
two pairs of real (unstable) eigenvalues occur in the linearized problem (|3|), while one pair of 
imaginary eigenvalues of negative Krein signature persists on the imaginary axis near A = ±i. 
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To compute the small negative eigenvalues of L_, we compute the leading-order potential 

otherwise and then proceed with the first-order potential: 



V r [ 0) = -1 for n = and K (0) 



h n = -l 

2, n = 

§, n = l 

-1, n = 2 

0, otherwise 



-1, n = -l 

-3, n = 1 

-2, n = 2 

0, otherwise 



Since the potential has no symmetry, we have to find the eigenvector of the eigenvalue problem 
© in the most general form if> n = Ae K<yn+2 \ n < —2, ip_i = B, ipo = C, ib\ = D, ip2 = E, 
and ^3 = Fe~ K ( n ~ 3 \ n > 3. Using the parametrization ([7]), we obtain at the leading order of 
0(6): 



Ee~ K , 



2)E = D, C + E=(e K + e- K -3)D, A = Be 



e K B = B + C. 



Since the equation at n = implies that C = —e(B + D) + 0(e 2 ), the chain of equations 
is uncoupled at the variables (D, E, F) and (A, B) at the leading order. Working with the 
chain for (D,E,F), we obtain the same cubic equation ([9]) with the same root e K = 2 + y/2, 
which gives the approximation [i = — e(l + + 0(e 2 ). Working with the chain for (A, B), 

we obtain the equation e K = 1, which hides a small root k = 0(e). To unveil this hidden 
eigenvalue, we have to extend the potential to the second order by 



t(2) 
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Using the parametrization (|7|) and the representation of the eigenvector in the form ip n 
Ae K ( n+3 \ n < -3, ^-2 = B, ^-i = C, ^ Q = D, Vi = ^, ^2 = i 7 , V's = G, and V>4 
we obtain at the leading order of O(e) + O(e ): 



H = Ge~ K , (e K -e)G = F, E+G = (e K + e~ K - 2 + ^ ) F. D+F 
and 



e K + e- R - 3 - y ) 1. . 



5e- K , ( e K --)B = C, B + D=(e K + e 



1 



3e 



where the approximation D = — e(C + -E) + 0(e 2 ) is sufficient for the purpose. Eliminating 
B, D, F, and G, we obtain two equations at the leading order O(e): 

-eE = (e K - 1 - e)C, -eC = {2e K + e" K - 3 - e)E 

Using the asymptotic expansion k = e«i + 0(e 2 ), we reduce the problem to a quadratic equa- 
tion for K\ with two roots K\ = 2 and Ki = 0. The non-zero root leads to the approximation 
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Figure 3: The same as Fig. [2j but for the discrete dark soliton with U+ 
[/_ = Z_ U {!}. 



.\{1}, U = {0}, and 



\x = — 4e 3 + 0(e 4 ). Each of the two roots above generates a negative eigenvalue 7 of L + L_ 
with the correspondence 7 = 2^(1 + 0(e)). Summarizing, the two pairs of real eigenvalues 
are given asymptotically by A = ±1.84776i/e and A = ±V 8e 3 , which are again found in Fig. 
[3] to be in very good agreement with the full numerical results. 

If U+ = Z+\{1}, Uq = {0, 1}, and U- = Z_, then Nq = o\ = 1 and one pair of real (unstable) 
eigenvalues occurs in the linearized problem ([3]), while two pairs of imaginary eigenvalues of 
negative Krein signature persist on the imaginary axis near A = ±i. To compute the small 
negative eigenvalue of L_, we use Wolfram's MATHEMATICA and compute the potentials 
of the eigenvalue problem up to the fourth order 



-1, n = {0,l} y{1) 
0, otherwise ' n 



-1, n = {-l,2} 
0, otherwise 
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n = {0,l}, 
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I' 



{-4,5} 
{-3,4}, 
{-2,3}, 
= {-1,2}, 
n = {0,l}, 
otherwise 



n = 
n - 
n - 

, n 



Using the parametrization (|7|) and the symmetry of the eigenvector ip n = ^_ n+ i, n G Z with 



e 



n ■ 
n 
n ■ 
n 
n 



1 

2, 
3, 
4, 
5. 



8 




0.02 0.04 0.06 0.08 0.1 % 

' x10" 



Figure 4: The same as Fig. EJ but for the discrete dark soliton with U + = Z + \{1}, Uq = {0,1}, 
and U- = Z_. 
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Figure 5: A typical example of the solution profile (top) and the spectral plane of its linearization 
spectrum (bottom) for the discrete dark soliton with U+ = Z + \{1, 2}, Uq = {0, 1, 2}, and U- = Z_ 
for e = 0.05. As predicted, the configuration is linearly stable for small e, bearing three pairs of 
imaginary eigenvalues (with negative Krein signature), but no real eigenvalue pairs. 



we obtain algebraic equations for coefficients A...F, which are solvable up to the the fourth 
order, subject to the characteristic equation e K = 1 + 2e 2 + 0(e 3 ). Therefore, k = 2e 2 + 0(e 3 ), 
such that // = — 4e 5 + 0(e 6 ). The small negative eigenvalue of L + L_ is thus approximated by 
7 = — 8e 5 + 0(e 6 ), while the pair of real eigenvalues of the stability problem (|3|) is given by 
A = ±\/8e 5 (l + 0(e)). The prediction for this small real eigenvalue, leading to a very weak 
instability in this case, is compared to numerical results in Fig. [U Once again, we report 
very good agreement between the two. 

• If U+ = Z+\{1,2}, Uq = {0, 1,2}, and U- = Z_, then iVo = and three pairs of imaginary 
eigenvalues of negative Krein signature persist on the imaginary axis near A = ±i. This is 
confirmed in Fig. \5\ showing a typical example of the discrete dark soliton and its linearization 
spectrum for e = 0.05. 
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In summary, we have offered a systematic way to assess the stability of discrete dark solitons 
in the prototypical dynamical lattice model of the DNLS equation. We have illustrated how the 
number of sign changes in the limiting configuration at the anti-continuum limit can be used to 
count the number Nq of small real eigenvalues of its linearization spectrum, when deviating from 
the anti-continuum limit. We have also associated the number of zeros in the limiting sequence 
with the number of isolated imaginary eigenvalues of negative Krein signature (which accounts for 
potential oscillatory instabilities for larger values of the coupling). In addition to this full charac- 
terization of the real and imaginary eigenvalues, we have offered a systematic approach towards 
computing asymptotic approximations of the relevant eigenvalues. In particular, we have connected 
small eigenvalues of operator L + L_ to the small eigenvalues of operator L_ and have developed 
perturbation series expansions in terms of the inter-site coupling constant. Within this method, 
relevant computations result in a finite-dimensional matrix problem. We have demonstrated this 
approach in a variety of configurations including the on-site and inter-site dark soliton structures 
of [31 [U [6], but also in multiple-hole configurations of p2], finding very good agreement between 
the analytical considerations and the full numerical results. It would be of particular interest to ex- 
tend relevant computations to higher dimensional settings, examining, for instance, the stability of 
discrete defocusing vortices in the two- or three-dimensional DNLS equations. Such considerations 
are deferred to future studies. 
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